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Abstract — We consider the optimal quantization of compressive 
sensing measurements following the work on generalization of 
relaxed belief propagation (BP) for arbitrary measurement chan- 
nels. Relaxed BP is an iterative reconstruction scheme inspired by 
message passing algorithms on bipartite graphs. Its asymptotic 
error performance can be accurately predicted and tracked 
through the state evolution formalism. We utilize these results 
to design mean-square optimal scalar quantizers for relaxed BP 
signal reconstruction and empirically demonstrate the superior 
error performance of the resulting quantizers. 

I. Introduction 

By exploiting signal sparsity and smart reconstruction 
schemes, compressive sensing (CS) 0], j2] can enable signal 
acquisition with fewer measurements than traditional sam- 
pling. In CS, an n-dimensional signal x is measured through 
m random linear measurements. Although the signal may be 
undersampled (m < n), it may be possible to recover x 
assuming some sparsity structure. 

So far, most of the CS literature has considered signal 
recovery directly from linear measurements. However, in many 
practical applications, measurements have to be discretized 
to a finite number of bits. The effect of such quantized 
measurements on the performance of the CS reconstruction 
has been studied in 0, §]. In 0-0 the authors adapt 
CS reconstruction algorithms to mitigate quantization effects. 
In JU, high-resolution functional scalar quantization theory 
was used to design quantizers for LASSO reconstruction |]9]. 

The contribution of this paper to the quantized CS problem 
is twofold: First, for quantized measurements, we propose 
reconstruction algorithms based on Gaussian approximations 
of belief propagation (BP). BP is a graphical model-based 
estimation algorithm widely used in machine learning and 
channel coding ifTOl . ifTTI that has also received significant 
recent attention in compressed sensing lfl2l . Although exact 
implementation of BP for dense measurement matrices is 
generally computationally difficult, Gaussian approximations 
of BP have been effective in a range of applications [ 13 1-[ 18 1. 
We consider a recently developed Gaussian-approximated BP 
algorithm, called relaxed belief propagation lfl6l . lfl9l . that ex- 
tends earlier methods 11151 . lfT8l to nonlinear output channels. 
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We show that the relaxed BP method is computationally simple 
and, with quantized measurements, provides significantly im- 
proved performance over traditional CS reconstruction based 
on convex relaxations. 

Our second contribution concerns the quantizer design. With 
linear reconstruction and mean-squared error distortion, the 
optimal quantizer simply minimizes the mean squared error 
(MSE) of the transform outputs. Thus, the quantizer can be op- 
timized independently of the reconstruction method. However, 
when the quantizer outputs are used as an input to a nonlinear 
estimation algorithm, minimizing the MSE between quantizer 
input and output is not necessarily equivalent to minimizing 
the MSE of the final reconstruction. To optimize the quantizer 
for the relaxed BP algorithm, we use the fact that the MSE 
under large random transforms can be predicted accurately 
from a set of simple state evolution (SE) equations [19|, |20|. 
Then, by modeling the quantizer as a part of the measurement 
channel, we use the SE formalism to optimize the quantizer 
to asymptotically minimize distortions after the reconstruction 
by relaxed BP. 

II. Background 

A. Compressive Sensing 

In a noiseless CS setting the signal x £ K™ is acquired via 
m < n linear measurements of the type 

z = Ax, (1) 

where A € R mx ™ j s the measurement matrix. The objective 
is to recover x from (z, A). Although the system of equations 
formed is underdetermined, the signal is still recoverable if 
some favorable assumptions about the structure of x and A 
are made. Generally, in CS the common assumption is that the 
signal is exactly or approximately sparse in some orthonormal 
basis i.e., there is a vector u = $!~ 1 x G R" with most of its 
elements equal or close to zero. Additionally, for certain guar- 
antees on the recoverability of the signal to hold, the matrix A 
must satisfy the restricted isometry property (RIP) [21]. Some 
families of random matrices, like appropriately-dimensioned 
matrices with i.i.d. Gaussian elements, have been demonstrated 
to satisfy the RIP with high probability. 

A common method for recovering the signal from the mea- 
surements is basis pursuit. This involves solving the following 



optimization problem: 

min ||\l/ _1 a;|| £i subject to z = Ax, (2) 

where || • is the ^i-norm of the signal. Although it is 
possible to solve basis pursuit in polynomial time by casting 
it as a linear program (LP) ll22l . its complexity has motivated 
researchers to look for even cheaper alternatives like numerous 
recently-proposed iterative methods 02, ED, 03. EH, ED- 
Moreover, in real applications CS reconstruction scheme must 
be able to mitigate imperfect measurements, due to noise or 
limited precision [3|, [5|, |6|. 

B. Scalar Quantization 

A quantizer is a function that discretizes its input by per- 
forming a mapping from a continuous set to some discrete set. 
More specifically, consider Appoint regular scalar quantizer 
Q, defined by its output levels C = {c»; i = 1, 2, . . . , N}, 
decision boundaries C K; i = 1, 2, ... , N}, and 

a mapping Cj = Q(s) when s e [6j_i,6j) [25]. Additionally 
define the inverse image of the output level c, under Q as a 
cell Q _1 (ci) = [6,_i,6j). For i = 1, if &o = — 00 we replace 
the closed interval [6 , 6i) by an open interval (b ,bi). 

Typically quantizers are optimized by selecting decision 
boundaries and output levels in order to minimize the dis- 
tortion between the random vector s 6 M. m and its quantized 
representation § = Q(s). For example, for a given vector s 
and the MSE distortion metric, optimization is performed by 
solving 

Q* =argmin£;{||s-Q(s)||, 2 2 }, (3) 

where minimization is done over all A-level regular scalar 
quantizers. One standard way of optimizing Q is via the Lloyd 
algorithm, which iteratively updates the decision boundaries 
and output levels by applying necessary conditions for quan- 
tizer optimality [25 1. 

However, for the CS framework finding the quantizer that 
minimizes MSE between s and i is not necessarily equivalent 
to minimizing MSE between the sparse vector x and its 
CS reconstruction from quantized measurements x |8|. This 
is due to the nonlinear effect added by any particular CS 
reconstruction function. Hence, instead of solving (|3), it is 
more interesting to solve 

Q* =argmmS{||x-x|| 2 }, (4) 
Q 1 J 

where minimization is performed over all A-level regular 
scalar quantizers and x is obtained through a CS reconstruction 
method like relaxed BP or AMP. This is the approach taken 
in this work. 

C. Relaxed Belief Propagation 

Consider the problem of estimating a random vector x £ 
R™ from noisy measurements y G M m , where the noise is 
described by a measurement channel p y \ z (y a | z a ), which acts 
identically on each measurement z a of the vector z obtained 
via (fl}. Moreover suppose that elements in the vector x are 



distributed i.i.d. according to p x (xi). Then we can construct 
the following conditional probability distribution over random 
vector x given the measurements y: 

Px|y | V) = -g Y[p x (Xi) Y[ Py\z (Va \ Z a ) , (5) 

i—1 a—1 

where Z is the normalization constant and z a = (Ax) a . By 
marginalizing this distribution it is possible to estimate each 
Xi. Although direct marginalization of p x \ y (x | y) is compu- 
tationally intractable in practice, we approximate marginals 
through BP 02, OH, 03- BP is an iterative algorithm 
commonly used for decoding of LDPC codes [11 1. We apply 
BP by constructing a bipartite factor graph G = (V, F, E) 
from (|5]l and passing the following messages along the edges 
E of the graph: 



P-ia fa) oc p x (a*) Y[ pUi ( Xi ) , (6) 




where oc means that the distribution is to be normalized so 
that it has unit integral and integration is over all the elements 
of x except a^. We refer to messages {pi^ a }(i,a)eE as vari- 
able updates and to messages {p a ^i}(i,a)eE as measurement 
updates. We initialize BP by setting p ^ a {xi) — p x (xi). 

Earlier works on BP reconstruction have shown that it 
is asymptotically MSE optimal under certain verifiable con- 
ditions. These conditions involve simple single-dimensional 
recursive equations called state evolution (SE), which predicts 
that BP is optimal when the corresponding SE admits a unique 
fixed point 031 . ll20l . Nonetheless, direct implementation of 
BP is still impractical due to the dense structure of A, which 
implies that the algorithm must compute the marginal of 
a high-dimensional distribution at each measurement node. 
However, as mentioned in Section Q] BP can be simplified 
through various Gaussian approximations, including the re- 
laxed BP method [15], [ 16 1 and approximate message passing 
(AMP) 03. H91 . Recent theoretical work and extensive 
numerical experiments have demonstrated that, in the case of 
certain large random measurement matrices, the error perfor- 
mance of both relaxed BP and AMP can also be accurately 
predicted by SE. Hence the optimal quantizers can be obtained 
in parallel for both of the methods, however in this paper we 
concentrate on design for relaxed BP, while keeping in mind 
that identical work can be done for AMP as well. 

Due to space limitations, in this paper we will limit our 
presentation of relaxed BP and SE equations to the setting in 
Figure [T] See [16| for more general and detailed analysis. 

III. Quantized Relaxed BP 

Consider the CS setting in Figure [T] where without loss of 
generality we assumed that ^ = I n . The vector x £ M. n is 
measured through the random matrix A to result in z g R m , 
which is further perturbed by some additive white Gaussian 
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Fig. 1: Compressive sensing set up with quantization of 
noisy measurements s. The vector z denotes noiseless random 
measurements. 



noise (AWGN). The resulting vector s can be written as 

s = z + r) = Ax + rj, (8) 

where {rj a } are i.i.d. random variables distributed as Af(0, a 2 ). 
These noisy measurements are then quantized by the A^-level 
scalar quantizer Q to give the CS measurements y € R m . 
The relaxed BP algorithm is used to estimate the signal x 
from the corrupted measurements y, given the matrix A, noise 
variance a 2 > 0, and the quantizer mapping Q. Note that 
under this model each quantized measurement y a indicates 
that s € Q~ 1 (ya), hence our measurement channel can be 
characterized as 



Py|z (Va I Z a ) 



iQ-Hva) 
for a = 1, 2, . . . , m and where 

1 



<j> (t - z a ; a 2 ) dt, (9) 



is Gaussian function 



S2irv GXP V 2v 



(10) 



Relaxed BP can be implemented by replacing probability 
densities in © and (0 by two scalar parameters each, which 
can be computed according to the following rules: 



~t + l _ 



p. 



J2b^a AbiUb^i 1 

v J2b^a ^-bi T b-yi Sb^a ^bi T b^i 



J2b^a ^bi T b^i J2b^a Ac 
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(12) 



b^a bi b—ti 



(13) 



= D 2 Va, A ajX ] A li^Ua + ° 2 I ' « 14i 

where a 2 is the variance of the components r\ a . Additionally, 
at each iteration we estimate the signal via 
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v^™ 42 _t ' v^' m 42 _t 



(15) 



for each i = 1, 2, . . . , n. 

We refer to messages {i'i-i-a, ^i^a}(i.a)e_E as variable up- 
dates and to messages {ii a _^, T _j.iW ) e .E as measurement 
updates. The algorithm is initialized by setting x^ a = x^t 



and t, 







of the prior p x (xi). The nonlinear functions F n and £; n are 
the conditional mean and variance 



F in (q,v) =E{x| q = q}, 
f in (q, v) = var {x | q = q} , 



(16) 
(17) 



where q = x + v, x ~ p x (a^), and v ~ Af (0, !/). Note that 
these functions admit closed-form expressions and can easily 
be evaluated for the given values of q and v. Similarly, the 
functions D\ and D% can be computed via 



Di (y, z, v) = - (z - F out (y, z, v)) 
n / - s _ 1 A ^out (y, z, v) 

D 2 [y, z, v) = - 1 

v V v 



(18) 
(19) 



where the functions F om and £ out are the conditional mean and 
variance 

F oal {y,z,iy) = E{z\zeQ- 1 (y)}, (20) 
E oat {y,'z,v) = var{z | z e Q^ 1 (y)} , (21) 

of the random variable z ~ J\f (z, v). These functions admit 
closed-form expressions in terms of erf (z) = f Q e~* dt. 

IV. State Evolution for Relaxed BP 

The equations (fTTT > — (fl~5T> are easy to implement, however 
they provide us no insight into the performance of the algo- 
rithm. The goal of SE equations is to describe the asymptotic 
behavior of relaxed BP under large measurement matrices. The 
SE for our setting in Figure Q] is given by the recursion 

vt+i = £m (fout (l3v t , a 2 )) , (22) 

where t > is the iteration number, (3 = n/m is a fixed 
number denoting the measurement ratio, and a 2 is the variance 
of the AWGN components which is also fixed. We initialize 
the recursion by setting vq = T; n i t , where fjnit is the variance 
of Xi according to the prior p x (xi). We define the function f n 
as 

£ m (v)=E{£ iB (q,v)}, (23) 

where the expectation is taken over the scalar random variable 
q = x + v, with x ~ p x (xi), and v ~ Af(0, v). Similarly, the 
function £ out is defined as 

fout {v, o- 2 ) = 7 - ( — \ : — jtti (24) 

E{D 2 (y, z, v + <7 z )\ 

where D 2 is given by ( [T9| > and the expectation is taken over 
p y \ z {y a I z a ) and (z, 2) ~ Af(0, P z {v)), with the covariance 
matrix 



P z (y) = 



^Tini, /3f init - v 
/3finit - V /3f init - V 



(25) 



finit where and fbit are the mean and variance 



One of the main results of [ 16 1, which we present below for 
completeness, was to demonstrate the convergence of the error 
performance of the relaxed BP algorithm to the SE equations 
under large sparse measurement matrices. Denote by d < rn 
the number of nonzero elements per column of A. In the large 
sparse limit analysis, first let n — > 00 with m — fin and 
keeping d fixed. This enables the local-tree properties of the 



factor graph G. Then let d —¥ oo, which will enable the use 
of a central limit theorem approximation. 

Theorem 1. Consider the relaxed BP algorithm under the 
large sparse limit model above with transform matrix A and 
index i satisfying the Assumption 1 of H16V for some fixed 
iteration number t. Then the error variances satisfy the limit 



R = 1 bits/component 



lim lim E ■ 

d— >oo n—. 



(26) 



'{l Xi x *L} 

where vt is the output of the SE equation \22\ . 

Proof: See HU. ■ 
Another important result regarding SE recursion in d22l is 
that it admits at least one fixed point. It has been showed that 
as t —> oo the recursion decreases monotonically to its largest 
fixed point and if the SE admits a unique fixed point, then 
relaxed BP is asymptotically mean-square optimal |[T6l . 

Although in practice measurement matrices are rarely 
sparse, simulations show that SE predicts well the behavior of 
relaxed BP. Moreover, recently more sophisticated techniques 
were used to demonstrate the convergence of approximate 
message passing algorithms to SE under large i.i.d. Gaussian 
matrices ifTSl , |fl9l 

V. Optimal Quantization 

We now return to the problem of designing MSE-optimal 
quantizers under relaxed BP presented in (|4). By modeling the 
quantizer as part of the channel and working out the resulting 
equations for relaxed BP and SE, we can make use of the 
convergence results to recast our optimization problem to 



Q SE = arg min \ lim v t \ 



(27) 



where minimization is done over all A^-level regular scalar 
quantizers. In practice, about 10 to 20 iterations are sufficient 
to reach the fixed point of D t . Then by applying TheoremQ] we 
know that the asymptotic performance of Q* will be identical 
to that of Q SE . It is important to note that the SE recursion 
behaves well under quantizer optimization. This is due to 
the fact that SE is independent of actual output levels and 
small changes in the quantizer boundaries result in only minor 
change in the recursion (see (BTTi). Although closed-form 
expressions for the derivatives of v t for large i's are difficult 
to obtain, we can approximate them by using finite difference 
methods. Finally, the recursion itself is fast to evaluate, which 
makes the scheme in (|27| | practically realizable under standard 
optimization methods like sequential quadratic programming 
(SQP). 

VI. Experimental Results 

We now present experimental validation for our results. 
Assume that the signal x is generated with i.i.d. elements from 
the Gauss-Bernoulli distribution 

N (0, 1 1 p) , with probability p; 
0, with probability 1 — p, 

where p is the sparsity ratio that represents the average fraction 
of nonzero components of x. In the following experiments we 



(28) 
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Fig. 2: Optimized quantizer boundaries for 1 bits/component 
of x. Optimal quantizer is found by optimizing quantizer 
boundaries for each j3 and then picking the result with smallest 
distortion 



assume p = 0.1. We form the measurement matrix A from 
i.i.d. Gaussian random variables, i.e., A a i ~ Af(0, 1/m); and 
we assume that AWGN with variance a 2 = 10~ 5 perturbs 
measurements before quantization. 

Now, we can formulate the SE equation (l22l and perform 
optimization (|27| ). We compare two CS-optimized quantizers: 
Uniform and Optimal. We fix boundary points bo — — oo 
and b jy = +oo, and compute the former quantizer through 
optimization of type (0). In particular, by applying the central 
limit theorem we approximate elements s a of s to be Gaussian 
and determine the Uniform quantizer by solving ©, but with 
an additional constraint of equally-spaced output levels. To 
determine Optimal quantizer, we perform (l27t by using a 
standard SQP optimization algorithm for nonlinear continuous 
optimization. 

Figure [2] presents an example of quantization boundaries. 
For the given bit rate R x over the components of the input 
vector x, we can express the rate over the measurements s 
as R s — /3R X , where /3 — n/m is the measurement ratio. 
To determine the optimal quantizer for the given rate R x 
we perform optimization for all /3s and return the quantizer 
with the least MSE. As we can see, in comparison with 
the uniform quantizer obtained by merely minimizing the 
distortion between the quantizer input and output, the one 
obtained via SE minimization is very different; in fact, it looks 
more concentrated around zero. This is due to the fact that by 
minimizing SE we are in fact searching for quantizers that 
asymptotically minimize the MSE of the relaxed BP recon- 
struction by taking into consideration the nonlinear effects 
due to the method. The trend of having more quantizer points 
near zero is opposite to the trend shown in [8 1 for quantizers 
optimized for LASSO reconstruction. 

Figure |3]presents a comparison of reconstruction distortions 
for our two quantizers and confirms the advantage of using 
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Fig. 3: Performance comparison of relaxed BP with other 
sparse estimation methods. 

quantizers optimized via (l22"l l. To obtain the results we vary 
the quantization rate from 1 to 2 bits per component of x, 
and for each quantization rate, we optimize quantizers using 
the methods discussed above. For comparison, the figure also 
plots the MSE performance for two other reconstruction meth- 
ods: linear MMSE estimation and the widely-used LASSO 
method (9|, both assuming a bounded uniform quantizer. The 
LASSO performance was predicted by state evolution equa- 
tions in |fl9l , with the thresholding parameter optimized by 
the iterative approach in [26|. It can be seen that the proposed 
relaxed BP algorithm offers dramatically better performance — 
more that 10 dB improvement at low rates. At higher rates, the 
gap is slightly smaller since relaxed BP performance saturates 
due to the AWGN at the quantizer input. Similarly we can see 
that the MSE of the quantizer optimized for the relaxed BP 
reconstruction is much smaller than the MSE of the standard 
one, with more than 4 dB difference for many rates. 

VII. Conclusions 

We present relaxed belief propagation as an efficient algo- 
rithm for compressive sensing reconstruction from the quan- 
tized measurements. We integrate ideas from recent general- 
ization of the algorithm for arbitrary measurement channels 
to design a method for determining optimal quantizers under 
relaxed BP reconstruction. Although computationally simpler, 
experimental results show that under quantized measurements 
relaxed BP offers significantly improved performance over tra- 
ditional reconstruction schemes. Additionally, performance of 
the algorithm is further improved by using the state evolution 
framework to optimize the quantizers. 
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